Time-dependent quantum many-body theory of identical bosons in a double well: 
Early time ballistic interferences of fragmented and number entangled states 
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A time-dependent multiconfigurational self-consistent field theory is presented to describe the 
many-body dynamics of a gas of identical bosonic atoms confined to an external trapping potential 
at zero temperature from first principles. A set of generalized evolution equations are developed, 
through the time-dependent variational principle, which account for the complete and self-consistent 
coupling between the expansion coefficients of each configuration and the underlying one-body 
wave functions within a restricted two state Fock space basis that includes the full effects of the 
condensate's mean field as well as atomic correlation. The resulting dynamical equations are a 
classical Hamiltonian system and, by construction, form a well-defined initial value problem. They 
are implemented in an efficient numerical algorithm. An example is presented, highlighting the 
generality of the theory, in which the ballistic expansion of a fragmented condensate ground state is 
compared to that of a macroscopic quantum superposition state, taken here to be a highly entangled 
number state, upon releasing the external trapping potential. Strikingly different many-body matter- 
wave dynamics emerge in each case, accentuating the role of both atomic correlation and mean-field 
effects in the two condensate states. 

PACS numbers: 03.75.Kk, 03.75.-b, 05.30.Jp, 03.75.Gg 



I. INTRODUCTION 

Understanding the quantum many-body structural 
and dynamical properties of the trapped gaseous Bose- 
Einstein condensate (BEC) lies at the heart of many 
experimental and theoretical research efforts worldwide. 
Remarkably sensitive matter-wave interferometry 0, 0, 
;3i], based upon an atomic BEC source, represents just 
one potentially useful experimental tool of technologi- 
cal import, while the demonstration of coherent macro- 
scopic superposition of millions of Bosc-Einstein con- 
densed atoms, may, one day, be turned from dream to 
reality, possibly answering questions fundamental to the 
theory of quantum mechanics Q. These examples, span- 
ning both practical and deeply fundamental extremes, 
have either been already realized [l], or the first fun- 
damental steps have been achieved [5| in the context of 
a BEC made of identical atoms confined to an external 
double-well trapping potential. In order to complement 
these experimental accomplishments, theoretical knowl- 
edge of the double-well condensate's structure and dy- 
namics, beyond that of simple models, is now in great 
demand. 

On the theoretical front, the time- independent quan- 
tum many-body structure of the double- well condensate 
is beginning to be explored with powerful first princi- 
ples methods adapted from the quantum-chemical theory 
of many-electron atomic and molecular systems [f| 0] at 
equilibrium and dissociation, reflecting the strong anal- 
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ogy between a rigorous description of molecular disso- 
ciation and that of BEC fragmentation. Such meth- 
ods involve the simultaneous variational optimization 
of both the Fock space expansion coefficients and the 
one-body wave functions (orbitals) underlying each Fock 
state Hil]. These methods build in, both, the full effects 
of the condensate's mean field and the correlations arising 
between atoms in different Fock states, and are capable 
of describing the many-body ground and excited states 
of the system at any barrier height for either symmet- 
ric or asymmetric trapping potentials. The unique role 
of spatial symmetry breaking in high-lying macroscopic 
self-trapped and superposition states of the double-well 
BEC has also recently been explored, for the first time, 
in detail, at a first principles level [lfj |. 

Building upon this detailed many-body structure 
comes the next frontier in the first principles theory of 
the double- well condensate: the dynamics. It is useful to 
learn from what has already been accomplished in the ex- 
plicitly time-dependent many-body approaches to atomic 
and molecular dynamics. Dynamical methods that are 
derivable from the principle of least action [ll|, E3] j the 
time-dependent variational principle [l3[ , or other equiv- 
alent variational approaches that provide a representa- 
tion of the time-dependent Schrodinger equation [l4j are 
regarded, by the authors, to be especially inspiring as 
they enjoy, by construction, various symmetries and as- 
sociated conservation laws, and lead to well-defined evo- 
lution equations that include the complete coupling of the 
dynamical variables within a chosen parameter space. In 
particular, description of the reactive chemical dynam- 
ics of general atomic or molecular systems, beyond the 
adiabatic and Born-Oppenheimer approximations, by a 
single determinant of coupled electronic parameters and 
nuclear coordinates and momenta [l5j represents one no- 
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table application of the time-dependent variational prin- 
ciple. Generalization of this work to the case of multiple 
determinants represents another. Also, recently, for- 
mally similar well-defined equations of motion describing 
the dynamics of the classical electromagnetic field cou- 
pled com plet ely to its atomic sources of charge and cur- 
rent (l7l. Il8l| have been developed with the aid of the 
time-dependent variational calculus. The Dirac-Frenkel 
variational principle [2(| applied to a multiconfigura- 
tional state vector for distinguishable bosons has been 
successful in characterizing, inter alia, the vibrational 
dynamics occurring in general polyatomic molecules [2l| . 
This same approach has been used to study the transition 
of a small and fixed number of indistinguishable bosons 
from a coherent to a fragmented and finally to a fermion- 
ized ground state by imaginary time propagation [22j . As 
noted in Ref . [lj| , the Dirac-Frenkel form of the varia- 
tional principle and the time-dependent variational prin- 
ciple used here in Sec. III. may yield superficially differ- 
ent looking evolution equations; however, if the variations 
are carried out fully and correctly, they will both give 
identical results. A notable time-dependent treatment of 
the double-well BEC, already in the literature, that intu- 
itively describes the dynamics of a double-well BEC by 
an improved two-mode model with time-dependent Fock 
space expansion coefficients but static Gross-Pitaevskii 
orbitals underlying each Fock state may be found in Ref. 




In the spirit of the above approaches to dynam- 
ics, a time-dependent multiconfigurational bosonic self- 
consistent field state vector is presented in Sec. II. that 
is capable of describing the many-body dynamics of a gas 
of identical bosonic atoms at zero temperature confined 
to an external trapping potential that can be continu- 
ously deformed from a single well to a double well. It 
is parametrized by complex-valued time-dependent Fock 
state expansion coefficients and atomic orbitals, which 
inherit the role of generalized coordinates in a nonlinear 
phase space. Equations of motion associated with this 
parametrized state vector are derived in Sec. III. by ap- 
plication of the time-dependent variational principle to 
the many-body Schrodinger Lagrangian. Properties of 
the evolution equations are discussed and a symplectic 
transformation is made to map the phase space coordi- 
nates to a real form that is amenable to numerical im- 
plementation. Lastly, in Sec. IV., numerical examples 
are presented in which the ballistic expansion dynam- 
ics of a fragmented ground state is compared to that of 
a macroscopic quantum superposition state (also called 
a Schrodinger cat state) of the BEC. These results are 
contrasted with a time-dependent Gross-Pitaevskii ap- 
proximation of ballistic BEC expansion. 



II. MULTICONFIGURATIONAL STATE 
VECTOR ANSATZ AND OVERLAP 

The most general time-dependent state vector repre- 
senting a gas of N identical bosonic atoms confined to an 
external double-well trapping potential at zero temper- 
ature should be flexible enough to describe the motion 
of atoms, in all possible arrangements, between the two 
trap minima, and, in addition, the dynamics of the state 
associated with each particular arrangement. Here it is 
assumed that two wells implies two Fock states. As the 
double- well experiments that are of primary interest to us 
involve symmetric trapping potentials [l], [2, 0, H| we fur- 
ther specialize to this symmetry (24[; therefore, within 
this restricted model space, a time-dependent state is 
sought that accounts for all possible arrangements of 
atoms. The collection of Fock states 

\$ Nl [<P]) = \NuN 2 ) = (b\)^(bl)^\vac}/^NjN^., 

(1) 

for Aq = 0, . . . , N and N 2 = N - Ni , with N k atoms re- 
stricted to two macroscopically occupied states, k = 1, 2, 
with the spatially orthogonal symmetric and antisym- 
metric one-body orbitals (/>fc(x, t) underlying each Fock 
state, provides such a basis. Therefore, in this restricted 
space, the most general completely symmetric many- 
body state may be written as the superposition 

JV 

|d(f); </>(*)) = £ \$ Nl [<j>(t)])d N dt), (2) 

N 1= 

where the vector d = d(i) = [do(t), . . . ,d,N(t)] and 
4> = 4>(t) = [<f>i (x, i), </>2(x, t)] are the collection of all 
N + 1 time-dependent expansion coefficients and sym- 
metric and antisymmetric space- and time-dependent or- 
bitals. Each term in the series, representing a particular 
arrangement of atoms between the two wells, is a called 
a configuration. This ansatz, which is a superposition of 
configurations, is a time-dependent generalization of the 
multiconfigurational state promoted in Refs. , where 
both the linear Fock state expansion coefficients and non- 
linear symmetrized product of orbitals within [</>]) 
are variationally optimized and mutually self-consistent; 
it is a time-dependent multiconfigurational bosonic self- 
consistent field (TDMCBSCF) state vector. Both the 
expansion coefficients d and the orbitals <j) of each config- 
uration |<E>jvi [4>]) are complex- valued and time dependent 
[25l | . It will be demonstrated that these parameters form 
a set of generalized coordinates whose canonically conju- 
gate momenta are proportional to their respective com- 
plex conjugates. Such a parametrized TDMCBSCF state 
spans a sufficiently general phase space to allow, upon 
application of the variational calculus, for a set of well- 
defined first order coupled nonlinear evolution equations 
with solutions enjoying a rich dynamics influenced by the 
full coupling between condensate's mean field and the 
quantum-mechanical correlations arising between atoms 
in different configurations. 
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A basic ingredient in the variational approach to deriv- 
ing dynamical equations is the overlap of two many-body 
states ([2]) . This multiconfigurational overlap involves the 
overlap of two single-configurational states \ &n ± [4>])- Pro- 
jection of |$jv"i[0]) onto the coordinate basis yields the 
single-permanental wave function 



$ Nl (l,...,N) = {l,...,N\$ Nl [<j>]) 
= 5{0i(l)».0i(JV 1 )0 2 (JVi + l) 



■MN)}, 



(3) 



where, for simplicity, the variable j = (x,-,t) for j — 

1, . . . , iV and S — (y/N\N 1 \N 2 l)- 1 J2p p is the s y m " 
metrizer with permutation operator P. The spatial over- 
lap of two such single-permanental wave functions is ex- 
pressible in a form that is reminiscent of that for deter- 
minants [26j, i.e., 

[</>*] 1$^ [</>]) 

= f{$ Nl [d>*]\l,...,N)(l,...,N\$ Nl [d>])dl---dN 
= • ■•fi(N l )<ft(N l + 1) • ■•<%(!*)} 

x 5{^i(l) • • • h(Ni)MNi + !)••■ MN)}dl ■■■dN 



t5{(0i|0i 



iVi!iV 2 ! 

X {(j>2 \<fa)Ni+l ■ ■ ■ (<H<fc)iv} 

= (iYi!7V2!) _1 permA(0*,0), 



(4) 



where a partial integration and S 2 = y^Nl/NjlN^lS has 
been used, and dj = d 3 Xj. The N x N orbital overlap 
matrix, denoted by A, reduces to 



A(0*,0) 



AjVixiVi 


(4>l < 




Ajvixjv 2 (^1! 




Aat 2 xJVi 




fe) 




h) 


AaTj x TVi 






OjVixJV 2 




OjV 2 xAfi 


Aat 2 xJV 2 (</ > 2) < 


h) 



(5) 



since the orbitals are orthogonal. The permanent of A is 
iViIA^!, which is related to the product of the dimensions 
of the N% x Ni upper left and N2 x N2 lower right blocks. 
Hence, the overlap of two permanents is the permanent 
of the overlaps and \^n x [<t>]) is normalized to unity. Note 
that time has not been integrated upon in Eq. so that 
the overlap remains time dependent. 

All multiconfigurational overlaps involve only single- 
configurational overlaps from the same configuration as 
different configurations are automatically orthogonal by 
symmetry. Together with Eqs. ©-(H]), the full multicon- 



figurational state overlap is given by 
5[d*,0*;d,0] = (d;0|d; < /,) 



N 

£ 

N 1= 



(^ 1 !Ar 2 !)" 1 ^ 1 [permA( < /»*,0)]d A r 1 . 

(6) 

Since the Fock state expansion coefficients satisfy 
X^tvi =0 I^JVi | 2 = 1, the time-dependent multiconfigura- 
tional state @ is unit normalized. 

In the following, derivatives of the overlap kernel S 
with respect to and <p* (and d and d*) will be 
needed to construct dynamical equations. Employing 
Laplace's expansion for permanents, i.e., perm A = 
^N\N[ [minor Aj^jy/ by expansion along column 
N[, these basic derivatives are 



(7) 



£ d *Ni N ^ldN l 
7Vi=0 



and 



d 2 S 
d<P* k 9<f>i 



N 

£ 

2Vi=0 
N 



(N 1 \N 2 l)~ 1 d* Ni 



d 2 [perm A] 



'■A'i 



(8) 



Ski)<l>k<f>*]dNi 

JVi=Q 



for k, I — I , 2. Other derivatives of import, such as 
d 2 S/dd*d(p or d 2 S/dd*dd, can either be obtained from 
these derivatives by complex conjugation or by differen- 
tiation with respect to d and d* . These latter derivatives 
are simple to compute from S. 



III. THE TIME-DEPENDENT VARIATIONAL 
PRINCIPLE 

Upon choosing a functional form for the time- 
dependent state vector, the time-dependent variational 
principle (TDVP) [l3j , which is derived from the princi- 
ple of least action applied to the many-body Schrodinger 
Lagrangian, generates a set of first order Hamiltonian 
equations of motion that form a well-defined initial value 
problem together with the initial values for the dynam- 
ical variables. The choice of parametrized state vector 
and its numerical representation are the only approxi- 
mations involved in this approach. Where the numerical 
representation of the state is systematically improvable, 
e.g., by enriching the computational basis with increased 
basis functions or grid points, the TDVP generates a hi- 
erarchy of systematically improvable equations of motion 
that form an approximation to the exact dynamical equa- 
tions. Indeed, in the limit of a complete basis, the TDVP 



4 



yields, in this case, the exact time-dependent many-body 
Schrodinger equation for N atoms in two macroscopi- 
cally occupied states. The TDVP equations and their 
solutions enjoy, by construction, a number of symme- 
tries and their corresponding conservation laws [13j |. and, 
in addition, preserve the orthogonality of the symmet- 
ric and antisymmetric orbitals <j>\ and </>2 without ad- 
ditional equations of constraint [24| , Furthermore, the 
TDVP equations exhibit the complete coupling between 
the dynamical variables allowed within the chosen state 
vector ansatz. This latter property, which is not nec- 
essarily guaranteed in other time-dependent approaches 
where intuition guides in the determination of couplings, 
is automatically achieved through the machinery of the 
calculus of variations; see, e.g., Ref. [27j |. 

The many-body Schrodinger Lagrangian L = 
L[d*, <j>*; d, <p] for a system of N interacting bosonic 
atoms, taken with respect to the multiconfigurational 
state @, is 

L = (d;<f>\ih(d/di) - H\d;<j)) 
= (^/2)(d;0| (d/dt) |d;0) - (ih/2){d;<f>\ (d/dt) |d;0) 



(d;#ff|d;</>), 



(9) 



where E = E[d* , </>*; d, </>] = (d; 0|.ff|d; 0) is the 
many-body energy and where the pure surface term 
— (ih/2)(d/dt)S has been added in the second line for 
symmctrization. Lagrangians differing only by a surface 
term always yield the same equations of motion and are 
called equivalent. Expansion of the ket total time deriva- 
tive (d/dt)— (<9/9£)£ and the bra total time derivative 
(d/dt)= (<9/<9£*)£*, with £ = [d,0], results in 



N 



- E 



dS u 



N 1= 



EE{ 

fc=l,2 q=l 



-<Pkq 



dS 



J kq 



>kq 



-E, 

(10) 



where the orbitals (j)k have been expanded onto an ar- 
bitrary real-valued basis {g q {x)}^ =1 of rank K so that 

0fc(x, t) ~ ^2^=1 9q( x ) ( t > kq(t). In this basis, the complex- 
valued expansion coefficients (j)k q = fikqit) are explicitly 
time dependent and inherit the role of dynamical vari- 
able. The details of this computational basis will be dis- 
cussed in greater detail in Sec. IV. 

The many-body Hamiltonian, which appears in the 



above Lagrangian, is given in second quantization by 

H = /#t(x,t)/i(x)#(x, t)d 3 x 

+ (l/2)/'l' t (x,t)* t (x',i)T/(x,x')*(x',t)*(x,t)d 3 xd 3 x' 
= ^ h kl (t)b\b l + (l/2) E V ki m n(t)b\b\ 



fcZ=l,2 



E 

klmn— 1,2 



(ii) 



and b\ 

creation operators making up the boson field operators 



where the b\ and b\ are basic boson annihilation and 



*(x,i) = 0i(x,i)6i+02(x,t)62 
*t( x ,t)=^(x,t)^ + ^(x,t)6+, 



(12) 



and h ki = ((j> k \h\<f>i) , and Vumn = (4>k4>i\V\4> m (l>n) are 
matrix elements of the one-body Hamiltonian /i(x) = 
(— h 2 /2m)'V 2 + T4 x t(x) and two-body atom-atom inter- 
action potential V(x, x']_= (ATrh 2 a/m)5(x — x') in the 
contact approximation [28j with repulsive s-wave scat- 
tering length a. In terms of this Hamiltonian, the energy 
is 



£[d*,0*;d,0] 



fci=l,2 



^fcZ7fci + (l/2) J3 VklmnXklnrm 
klmn— 1,2 

(13) 

where j^i and Tkimn are Fock space matrix elements of 
the parametrized complex-valued one- and two-body re- 
duced density matrices 7(1,1') and T(l, 1'; 2, 2'). Their 
functional form will be specified in Sec. III. C. Note that 
the many-body energy expectation (|13[) , which is a func- 
tion of the dynamical variables, is multiconfigurational; it 
is not the energy of a single configuration but, rather, rep- 
resents the energy of all interacting configurations within 
the restricted Fock space. 



A. The variational calculus 

Carrying out the full calculus of variations or the 
principle of least action [ll|, [l^] on the many-body 
Schrodinger Lagrangian (jTUJ) generates a set of well- 
defined equations of motion that approximate the exact 
time-dependent many-body Schrodinger equation for a 
system of N identical bosonic atoms at zero tempera- 
ture in an external trapping potential V eX f Using the col- 



lective notation £ = [d, <j> 
these TDVP equations [11 



for simplicity and generality, 
I12JI are now constructed from 



the Lagrangian L = (ih/2)[(dS/d£)£ - (8S/dC)C 
in Eq. (SJJ. 

The momentum conjugate to £ is 

so that its total time derivative becomes 



r)T 

(d/dt)^ = (ih/2) 
= (ih/2) 



— £ H £ 



[8 2 S 



d 2 S 



d£dC d£*<9£ 



as 

5£ 
t 



E 



(14) 



(15) 



5 



Lastly, the derivative of L with respect to £ itself is 



f = («/2) 



d 2 S 



a 2 5 



-r 



W (16) 



The TDVP equations, which are of first order in time, 
may then be built up from (d/dt)8L/8£ = dL/d£ and its 
complex conjugate, i.e., 



.. d 2 S ■ dE 



d?ac 



i 



where it has been assumed that the mixed partial deriva- 
tives of 5 commute with respect to £ and £*. It is these 
derivatives of the overlap that were computed, in part, in 
Sec. II. Eqs. (jTTJ) form a classical Hamiltonian system. 

Having given a general derivation of the Hamilton's 
(TDVP) equations for the collective dynamical variables 
£ and £*, the complete set of TDMCBSCF dynamical 
equations are now presented. Following the same ma- 
chinery as previously described, they are: 



N 



d 2 S 



K 



7^ dd Nl dd* 
i 1 

N o2 



4i + EE 



d 2 s 



dd* Ni dd; 

8 2 S 
dd* 



N{=0 Nl 



' 1 ' 

<*[£ s ' s 



N , =0 -^ kq oa N 

N ^2 



tq dd N 



K 



-d N , + E E 



4 ; + EE 

1=1,2 q' = l 
K 

: d N , + e E 

i=l,2 q' = l 



d 2 S 




^Iq' 


dE 


dd^dcj: 


W 




dd Nl 


d 2 S 




kq> 


dE 


dd* Ni d4 


— < 

lq> 


8d* Nl 


d 2 S 






8E 


d<j) kq d(i 




4>lq' 


8(f>kq 


d 2 S 




<t>lq> 


8E 


dj>t q d 4 







(18a) 
(18b) 
(18c) 
(18d) 



r 



These equations are closed and coupled, as is already 
evident, and will be shown to be highly nonlinear. Once 
d, d*, </>, and <j>* are specified at some point in time, then 
Eqs. (|18[) form a well-defined initial value problem. The 
resulting dynamics unfolds in a generalized phase space 
of dimension 2(N + 1) + 4K whose coordinates are the 
generalized positions and momenta d, d*, <j>, and <fi* . 

By restricting the time-dependent state vector to 
a single configuration, the TDVP would generate a set 
time-dependent SCF or Hartree-Fock equations for iden- 
tical bosons in two orbitals. Such a derivation has already 
been performed in the identical fermionic case for arbi- 
trarily many orbitals [1^, [3(| ■ An analogous set of time- 
dependent mean-field equations has also recently been 
derived for identical bosons in arbitrarily many orbitals 



B. Equations of motion in canonical form 

The TDVP equations of motion presented in Eqs. (TT5)) . 
which contain the full coupling allowed within the ansatz 
([2]) between the Fock state expansion coefficients d of 
each configuration |3>iVi [</>]) an d the underlying orbitals 
</>, may be collected into matrix form to simplify and, 
simultaneously, highlight their structure. Introducing the 



notation 



C, 



d 2 S 



(19) 



Eqs. HH]) become: 



- ihC dd d* - ihC d(j} <p* = dE/dd (20a) 

ihC Ad d + ihC d 4,4> = dE/dd* (20b) 

-ihC^d* - ihC^<p* = 8E/84> (20c) 

ihC^d + ihC^ = dE/dcf)* (20d) 

These TDMCBSCF equations may be organized into ma- 
trices according to 



Cdd 





~ °dd 


Cd<t> 






°d</> 





^d^ 









" d " 




' dE/dd* ' 


d* 




dE/dd 







dE/dcf)* 


0* 




_ dE/dc/} _ 



They are of the canonical form 

dH 



ujrj 



dr] 



(21) 
(22) 
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where 77 is a column vector whose entries are the dynam- 
ical variables, i.e., 77 = [d, d*, <fi, cf>*}. The matrix 



D J 
Jt M 





Cdd 





Cd0 










_ °dd 





™°d0 


ih 




















°d0 





C00 










d</> 





/I* 

-°00 J 



(23) 

which multiplies the column of velocities on the left hand 
side of Eq. ([2H)l , is called the symplectic form [ll|, [12] ■ It 
is, in this case, a nonlinear complex- valued function of the 
dynamical variables, i.e., lu = w[d*, </>*; d, 0], and defines 
the symplectic structure of the generalized phase space of 
the dynamical system [32[ • It contains the unique matrix 
elements: 



[Cdd] NiN{ 
[Cd(j,]N 1 k 



Cd Nl <t> k = N k <j)* k d Nl 



(24a) 
(24b) 
(24c) 



JV 



E d ki N i s > 

Ni=0 



Nt(N k 



dE 



For later simplification, the Fock state expansion coeffi- 
cient sector of w is blocked into the 2(N + 1) x 2(iV + 1)- 
dimensional matrix D, the upper right rectangular 2(N+ 
1) x 4i4T-dimcnsional matrix J provides the coupling be- 
tween the expansion coefficients d and the underlying 
orbitals <fi, and the orbital sector of ui is blocked into the 
AK x 4i^-dimensional matrix M. It is noted, that much 
of the nontrivial coupling between dynamical variables is 
contained in u>. The remaining coupling occurs in forces, 
as will now be demonstrated. 



Generalized forces 



Appearing on the right hand side of Eq. (fT8|) are the 
generalized forces 



dd* Ni (t) 



= \Wi + 2)(JVi + 1){N 2 - l)N 2 (l/2)V 2211 (t)d Nl+2 {t) 



+ {V(Ni + l)N 2 [h 21 (t) + 7ai U (t)JVi + V 2221 (t)(N 2 - l)]}d Nl+1 (t) 

+ {feii(t)JVa + h 22 (t)N 2 + (l/2)ry m2 (t) + VtiM^NtNi + (l/2)V lin (t)(Nf - N,) + (l/2)V 2222 (t)(iVf - N 2 )}d Nl (t) 
+ {y/Ni(N 3 + l)[h 12 (t) + V ni 2(t)(N 1 - 1) + Vu^N^dN^t) 
+ VMNi - 1)(N 2 + 1)(N 2 + 2)(l/2)V 1122 (t)d Nl - 2 (t) 

(25) 



dE 



[h(x)j kk (t) + V fefc (x, t)T kkkk (t) + Vfc/fc/(x, t)T kk , kk ,(t) +V kk ,{x,t)T kkkk ,(t) + Vfc<fe(x, f)IV fefefe (f)]^ fe (x, i) 

+ [/i(x)7fc fc / (t) + V fcfc (x, t)T kkkk > (t) + Vk>k'(x-,t)Tkk'k>k>(t) + Vfefe'(x,t)r fc fefe'/b'(*) + Vfe/fc(x,t)r fc fc/j/fc(t)]^(x,t) 

(26) 



r 



for iVi = 0, . . . ,iV and fc ^ /V = 1, 2. Here, the com- 
plex conjugate forces are omitted for brevity. The forces 
are written in terms of the Fock space matrix elements 
7^; and T k i mn of the complex-valued one- and two-body 
reduced density matrices 7(1,1') = (# t (l)*(l / )) and 
r(l,l';2,2') = (*t(i)#t(i')#(2)#(2'j) ll, i| taken 



with respect to the TDMCBSCF state ©. That is 
7(1,1') = {d, <&[</>*] I* 1 (l)tf(l') \d,m) 



N 



fel=l,2 N!N[=0 

= E ^(1)7^^(1') 



fcZ=l,2 



(27) 
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r(l,l';2,2') = (d,$[01|* t (l)v|/t(l')*(2)v|/(2')M,$[0]> 



E E rfk< $ iv I [0*]ist^s TO s„i$ JVi [0])d iV ^ m (2)^(2') 

klrnn 



(28) 



r 



for fc, l,m,n= 1, 2. Note that both 7 and T are nonlinear 
functions of the dynamical variables. The bosonic field 
operators $(x, t) and ^(x, i) appearing in these equa- 
tions, satisfy standard boson commutation relations. 

One-body matrix elements of the atom-atom 
interaction potential are given by Vfc;(x, t) = 
J 4>* k (x.',t)V(x,x.')(j>i(x',t)cl 3 x'. The symmetrized, direct 
plus exchange, matrix elements of V(x, x') are written 
as Vkimn = Vkimn + Vkinm for simplicity. All of the above 
matrix elements, both in Fock space and in configuration 
space, enjoy symmetry relations among the indices, e.g., 



v kl = vr k , v kl 



Vlkr, 



V* = V* 
mnkl nrnlk' 



and jkk' — Ik'k (fc 7^ fc' = 1,2). The symmetries of the 
Fock space matrix elements of the Hermitian two-body 
reduced density matrix 



r = 



^kkkk ^kkkk' 




k 


*-kkk'k' 


^kk'kk' 




k 


r kk' k'k' 




^k'kk 


k 


Fk'kk'k' 








^k'k'k'k 



(29) 



reveal that only six out of the ten matrix elements in 
the upper triangle of Eq. (f2l)|) are independent; they are 

Tfefefefe, Tk'k'k'k'-, ^kk'kk', ^kkk'k', ^kkkk' , and Tkk'k'k'- A 

similar situation exists for fermions with the labels k and 
kl replaced by holes and particles; see, e.g., Ref. [if] ]. 

It is now evident that the TDVP approximation of the 
exact time-dependent many-boson Schrodinger equation 
leads to dynamical equations that are highly nonlinear 
both through the symplectic form lj and through the 
forces dH/drj, which are themselves dynamical functions. 
Linearization of these multiconfigurational TDVP equa- 
tions forms the basis for a bosonic multiconfigurational 
random phase approximation; see, e.g., Ref. [34j | for the 
multiconfigurational fermionic case. 



D. Symplectic transformation to real- valued 
dynamical coordinates for numerical efficiency 

In this work, preference is given to numerical rou- 
tines that involve real-valued equations and their solu- 
tions rather than their complex- valued analogs. In effort 
to achieve this partiality, a symplectic transformation is 
constructed to map Eqs. (|23[) into manifestly real form. 
A symplectic (canonical) transformation is a mapping 



from a set of dynamical variables to a new set of dy- 
namical coordinates that is constructed in such a way as 
to preserve the symplectic structure of the generalized 
phase space (32j |. The transformation 



T = 



I 

71 



X 
Y 



(30) 



where both the 2(N + 1) x 2(N + l)-dimensional X and 
the 4K x 4A"-dimensional Y matrices are of the block 
form 



diag(l) diag(l) 
diag(-i) diag(i) 



(31) 



is constructed to map (£,£*) — * (R- e {£}: I m {£}) an dj si- 
multaneously, preserve the structure of the dynamical 
equations. That is, under T, the equations of motion 
(|2"3")l become 

~* TujT ^ = T ^- (32) 
The latter equation may be written more compactly as 

(33) 



(TujT^)f) = lot) = 



where fj = Tr\. These transformed dynamical equations 
are manifestly real valued. The new phase space coordi- 
nates and derivatives are 

fj = [d^d 1 ,^,^] = \/2[Re{d},Im{d},Re{(/)},Im{0}] 

(34) 

and 



d/dfj = [d/dd R , d/dd 1 , d/dcf> R , d/dcfr 1 } 



(35) 



and the new antisymmetric (or skew-symmetric) sym- 
plectic form is 



U= 2 



xdx^ Ajyt 



(36) 



A detailed illustration of each block of Cj is now provided. 
The upper left 2{N + 1) X 2(iV + l)-dimensional block 
becomes 



1 



XDX^ 



diag(-l) 
diag(l) 



(37) 
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while the upper right rectangular 2(N 
dimensional block is 



1) x 4S- 



and the lower right 4K x 4S-dimensional block is given 
by 



-[XJY% i;lq , 



dm 



Re I 



-N,Tm(d Nl 4>i q >) -JVjRefdjv-!^) 
Ni~Re{d Nl (f>* q ,) -Nilm(d Nl ^ ql ) 



1 



YMY^ 



A B 
-B T A 



(39) 



in terms of the 2K x 2S-dimensional sub-blocks 



(38) 



A kn , 



kq;lq' 



-T nn Im((f>i g (f>* qf ) -T 1212 lm(<p lq (j)^ q ,) 

-T 2121 lm{(j) 2q (j)* lql ) -r 2 222lm(0 2g ^2 g ') 



N 



En 1= o \d Nl \ 2 N2N 1 lm( ( j )2q r iql ) - En 1= o \dm\ 2 N 2 (N 2 - l)Im(<fc^v) 



N 



Bl,„: 



kq\lq' 



[7n<W +r m iRe(^iq<?!>^,)] 

-r 2 mRe(<K^) 



-ri21 2 Re(01 g 02g') 
[722<W +r 2 222Re(</>2 g '/>2 9 ')] 



E£ 1= o \dN l ?[Ni8 qql +N 1 {N 1 - l)Re{Ml q ,)] 



E£ 1= o \d Nl \ 2 [N 2 5 qq > + N 2 (N 2 - l)Re(0 2 ^,)] . 



(40) 



In both Eqs. ([38)1 and (|39)) . the contraction is on I = 1, 2 
and = 1, . . . , K. The lower left rectangular AK x2(N+ 
l)-dimensional block of Q is the negative transpose of Eq. 



The transformed real-valued forces are simply related 
to the real and imaginary parts of Eqs. ([23)) and (|26p. 
That is 



dE 
dd* 
dE 
3d 7 
dE 
d^ 
dE 
d^ 1 



1 
1 

1 



as dE 

dd* + ~dd 

. as 

dE 



dE 
~dd 



V2ldcj>* 
1 

vi 



as 
.as .as 



(41) 



As before, once the initial values for 77 arc specified, then 
the real-valued TDMCBSCF equations ((33]) form a well- 
defined and numerically efficient initial value problem, 
that is equivalent to the complex- valued equations (|2ip . 
and contain the complete coupling between the dynami- 
cal coordinates parametrizing the TDVP state vector (J2J) . 
The real- valued phase space, in which the solution of Eq. 
([33]) evolves, is endowed with the Poisson bracket 



as ^dG 
dfj 



(42) 



many-body Schrodinger equation takes on the classical 
Hamiltonian form 



V = {V,H}- 



(43) 



These dynamical equations have been completely imple- 
mented in an extensive computer program from which 
the first numerical results will now be presented. 



IV. NUMERICAL IMPLEMENTATION 

Together with an initial value for the vector 77, the 
above set of first order Hamiltonian (TDVP) dynamical 
equations is well defined. A simplistic quadrature scheme 
such as 

as 

rj(t + e)= rj{t) + efj{t) = n(t) + eu-\t)—— (44) 



dv(t) 



may be written to integrate the state vector forward in 
time from t to t + e 35]. This does require that the sym- 
plectic form u> be inverted at each time step; however, 
it has been found that, in practice, ui is readily inverted 
with standard AX = B inversion routines [3g]. These 
routines are easily incorporated into numerical integra- 
tors that are more sophisticated than Eq. ([44)) . This 
matter will be discussed in greater detail in the following. 



Preliminaries 



so that, together with the symplectic structure generated 
by {•, •}, the TDMCBSCF representation of the exact 



The TDMCBSCF evolution equations ([55)1. generated 
from the TDVP, have, so far, been represented in an ar- 



9 



bitrary basis of rank K of functions g q {x). In practice, 
these equations are solved using a fast Fourier transform 
based pseudospectral grid method [H, Hil in quasi-one 
dimension (iol. l4l|. The Fourier grid basis is constructed 
so that the basis functions g q (x) enforce the appropriate 
boundary conditions on the orbitals at the grid bound- 
aries. For the box boundary conditions employed here, 
these basis functions are the Fourier sine functions. It 
is found that the solutions of Eq. (j33|) are adequately 
converged with K — 2 9 fixed grid points in quasi-one di- 
mension. For the short time dynamics presented in the 
following, this value of K is sufficiently large enough to 
avoid the case where the BEC density is appreciably dif- 
ferent from zero at the box edges. Propagation to longer 
times, on the order of 50 ms or larger, may require the 
use of more grid points. The extension to three dimen- 
sions poses a challenge that may be overcome, without 
difficulty, with increased computing time. 

The equations of motion (I33p are integrated forward 
in time in real space, where a Fourier transform to k- 
space is made at each time step to diagonalize the spatial 
derivatives appearing in the generalized forces (|25p and 
(f2l)|). These terms are then transformed back to direct 
space before time propagation. Inversion of the symplcc- 
tic form Co is achieved, in direct space, with standard lin- 
ear algebra routines [36| and is incorporated into a vari- 
able step size fourth order Runge-Kutta method adopted 
from Ref. [33] - Difficulties arising from the inversion of 
Co, which is on the order of several thousand by several 
thousand, have yet to be encountered; however, analyti- 
cal and numerical methods do exist to correct this situa- 
tion should it arise in the future [lU, [52] • These methods 
rely on Darboux's theorem (see, e.g., Ref. [32}) to find a 
local coordinate chart in a curved phase space where Co 
takes on the constant form 



In a flat phase space, it is always possible to find a global 
system of coordinates where this result holds. Since the 
TDMCBSCF phase space is curved, transformation of Co 
to constant form is only possible locally. 

Values of the physical parameters such as the atomic 
mass m, scattering length a, and oscillator length f3 — 

h/mco, which enter the Hamiltonian are taken 

from Ref. H, as is the functional form of the external 
trapping potential V ey ±. These values are scaled, in the 
quasi-one-dimensional approximation (40l . l4lj , so that the 
product aN, with N — 100, is consistent with the 23 Na 
double-well interference experiments [3, 0, El] performed 
at MIT. It is these experiments, and those in Heidelberg 
Q, that have inspired the following numerical examples. 
Note that the TDMCBSCF theory is not limited to the 
quasi-one-dimensional approximation, the contact inter- 
action approximation, nor to the particular choice of V ox t 
that is used in the following. 



B. Numerical Results 

Imaging of typical laboratory BECs is commonly per- 
formed by releasing the condensate's trapping potential, 
and, subsequently, allowing the degenerate gas to ballis- 
tically expand to a point where it is large enough to be 
photographed. When atoms are, e.g., condensed into a 
double-well trap, release of the trap enables the initially 
left- and right-localized atomic clouds, which may or may 
not be coherent, to expand and eventually overlap due to 
their mutually repulsive interactions [43j. The details of 
the overlap dynamics depend strongly upon the initial 
phase coherence of the BEC and, additionally, upon the 
degree of atomic correlation between atoms in different 
Fock states. 

Using the TDMCBSCF theory, we present the ballis- 
tic expansion dynamics of two initially stationary states 
of the double-well condensate of N identical atoms at 
zero temperature, where, in each case, the left- and right- 
localized condensate moieties are: 

(1.) initially stationary and, consequently, phase coher- 
ent across both wells 

(2.) initially phase offset in each orbital by (tt/A)/N\ in 
the left hand well 

The two stationary states of interest are the fragmented 
ground state and a high-lying macroscopic quantum su- 
perposition state of the BEC. The latter is a highly en- 
tangled number state and is also called a Schrodinger cat. 
Both are obtained, initially, from full time-independent 
MCBSCF calculations; see Refs. (EGS- Following the 
release of the trap at t = 0, neither of these states are 
stationary. It is their many-body density 

p(x, t) = 7(x, t; x' = x, t = f) 

N 

= E \d Nl (t)\ 2 [Ni\M^t)\ 2 +N 2 \M^t)\ 2 } 

Ni=0 
N 

+ J2 \M(^2 + i)d^ 1 (t)0;(x,t)^ 2 (x,t)d J v 1 _i(t) 

N-l 

+ E y/W- + l)N2d* Nl (t)r 2 (^t)M^t)dN 1+ i(t), 

Ni=0 

(46) 

orbitals (x, t) , and expansion coefficients (t) that 
are followed as functions of space and time. 

1. Ballistic expansion dynamics of an initially fragmented 

state 

The fragmented TDMCBSCF state is initially single- 
configurational, i.e., it involves only the single configura- 
tion 

\d;<f>) = \$ N/2 [<p]) = \N/2,N/2) (47) 
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FIG. 1: (Color online) Ballistic expansion dynamics of the initially fragmented ground state density of the double-well conden- 
sate computed within the TDMCBSCF theory. The upper three panels display the dynamics of an initially stationary state of 
the BEC, while the lower three panels display that of an initially (n/4) /(N/2) orbital phase offset between the left and right 
wells. Time increases from left to right, with snapshots taken every 0.25 ms, so that a history of the dynamics is captured at 
1 ms and 2 ms after release of the trap at t = 0. It might appear that only three snapshots exist in frame one, but the color 
version will show that the density at t = 0.25 ms overlaps the initial density. Colors alternate to aid in visualization. Here, the 
short time dynamics reveals clearly visible interference between the left and right pieces of an initially fragmented condensate. 
The location of the interference fringes are not random since, in TDMCBSCF theory, the initial fragments have a well-defined 
phase relationship. The insert in panel three is taken, with permission pending, from the early BEC interference experiments 
performed at MIT where two independent condensates were allowed to ballistically expand and overlap; see Ref. [43]. Note 
that the experimentally observed interference fringes do not have 100% visibility in agreement with our TDMCBSCF theory. 
It will be demonstrated that Gross-Pitaevskii theory vastly overemphasizes the visibility of these fringes. 



with N/2 atoms in each well for N = 100. No corre- 
lations exist at t = between the atoms, although a 
specific phase relation must be chosen between left and 
right fragments; correlations do arise as the dynamics un- 
folds. Fig. [1] displays the short time expansion dynamics 
of the fragmented ground state density. The dynamics 
of an initially stationary state is presented in the upper 
three panels, while that of a condensate having an ini- 
tially (-k/A)/Ni orbital phase offset between the left and 
right wells is displayed in the lower three panels, where, 
in this case N\ = N/2. Time increases from left to right 
with each snapshot separated by 0.25 ms, so that the 
fragmented ground state dynamics is captured at 1 ms 
(panel one) and 2 ms (panel two) after the release of the 
trap at t = 0. For clarity, the third panel redisplays the 
density snapshot at 2 ms in the foreground and the ini- 
tial state density in the background. In each panel, the 
density at the final integration step is plotted in black. It 
can be seen that an initially fragmented BEC state does 



produce clearly visible interference fringes upon overlap 
of the left and right fragments. In an experiment with 
random relative phase between two independently pre- 
pared condensates, the location of the interference fringes 
would be random 44, 45, 46]; however, in this theoretical 
model, there is always a well-defined phase relationship 
between the two moieties, and, consequently, the inter- 
ference fringes are not randomly located. Note that the 
specific but natural choice of zero phase offset between 
left and right fragments has been assumed in the upper 
three panels of Figs. [1] and [51 The insert in panel three 
of Fig. [T]is taken from the early BEC interference exper- 
iments performed at MIT where two initially indepen- 
dently prepared condensates were allowed to ballistically 
expand, overlap, and interfere; see Ref. [43|]. Note that 
the visibility of the experimentally observed interference 
fringes is not 100%, which is in agreement with our fully 
variational TDMCBSCF theory; a similar result is found 
from a first order perturbative treatment in Ref. [J?) . In 
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the following, it will be shown that Gross-Pitaevskii the- 
ory, which is a single-configurational approach, greatly 
overemphasizes the visibility of these interference finges. 

The underlying orbitals 0&(x, t) and expansions coef- 
ficients djvi(0) corresponding to the initially fragmented 
density in Fig. [TJ were obtained, from first principles, 
with our time- independent MCBSCF theory [f|. In gen- 
eral, the construction of an arbitrary self-consistent ini- 
tial state can be difficult. Here, we rely on the fact 
that, within our restricted model space, the state of a 
left and right fragmented BEC can be re-expanded as 
a complicated (but equivalent) multiconfigurational su- 
perposition of symmetric and antisymmetric states with 
varying numbers of atoms in each. Since the TDMCB- 
SCF state is an arbitrary vector in the space spanned 
by either the left /right or symmetric / antisymmetric rep- 
resentation, it may be equivalently represented in cither 
basis. This fact is used to ensure that the underlying 
symmetric and antisymmetric orbitals c^fc(x, t) used here 
actually correspond to an initially fragmented state. 

Fig. [2] presents a series of snapshots of the orbitals 
<^fc(x, t) and expansion coefficients cZ/v^t) of the initially 
fragmented TDMCBSCF state |d;0) = \N/2,N/2). 
Each snapshot is separated by 0.25 ms. The expansion 
coefficients, which represent the distribution of atoms be- 
tween the left and right wells, initially correspond to a 
sharp \N/2,N/2) single-configurational state. Following 
release of the external trapping potential, it can be seen 
that d,N 1 (t) broadens in time and begins to approach a 
binomial distribution. It is well known that, in two-state 
multiconfigurational models, the expansion coefficients 
corresponding to the coherent symmetric ground state 
take the form of a binomial distribution peaked around 
the \N/2, N/2) configuration; see, e.g., Ref. [H[ and Ref. 
Q. In light of this fact, we note that the state of the 
BEC is moving from a single-configurational description 
that is initially left and right localized in each well to 
a multiconfigurational delocalized description. In an ap- 
propriate basis of symmetric and antisymmetric states, 
the latter delocalized multiconfigurational state reduces 
to the single-configurational Gross-Pitaevskii (GP) solu- 
tion. The initial fragmented state is well described by 
Hartree-Fock (HF) theory with left and right localized 
orbitals. Similarly, the coherent delocalized ground state 
is well described by GP theory with a single symmetric 
orbital. However, it is impossible, within either of these 
single-configurational approaches, to transition from the 
first to the second description. For example, by allowing 
a left localized and a right localized orbital to expand 
and eventually overlap in time-dependent HF theory, no 
interference occurs at any time since HF theory is a non- 
interacting single-particle theory. Furthermore, the final 
state of the system, following overlap, does not settle 
down into a symmetric coherent state solution of the GP 
equation. It is especially important to note that a full 
time-dependent multiconfigurational theory, such as the 
TDMCBSCF approach, is necessary to correctly describe 
the interaction physics and properly transition between 



the two single-configurational pictures [49| ■ 

2. Ballistic expansion dynamics of an initially number 
entangled state 

In order to explore the effects of initial atomic correla- 
tion on the dynamics of the double-well condensate, re- 
sults are presented for the ballistic expansion of a macro- 
scopic quantum superposition state, which is double- 
configurational. It represents the simplest possible ex- 
ample where correlations are present at t = 0. To this 
end, we choose a (superposed) number entangled state 
vector of the form 

|d; 0) = d Nl \$ Nl [0]) + d N2 |$at 2 [0]) 

= (i/^)|JVi,jv 2 ) + (i/V2)|jv 2j JVi) i4s) 

with Ni = 90 and N 2 = 10 for N = 100 atoms. Fig. 
121 displays the short time expansion dynamics of the su- 
perposition in Eq. (|48|) after the trapping potential is 
released at t = 0. The dynamics of a superposition that 
is initially stationary is presented in the upper three pan- 
els, while that of a superposition with underlying orbitals 
that are initially phase offset by (tt/4)/(Ni — N2) between 
the left and right wells is presented in the lower three 
panels, where, in this case iVi = 90 « N 3> N2 and, 
consequently, we take {ir/&)/(N\ - N 2 ) « (tt/4)/JVi w 
(ir/4)/N. Therefore, the orbital phase offset of (ir/4)/N 
leads to a macroscopic phase offset of 7r/4 in Eq. (|48p . 
Snapshots are taken, from left to right, every 0.25 ms so 
that the dynamics is captured at 1 ms, 2 ms, and 3 ms 
following the trap's release. For clarity, the history end- 
ing at 3 ms shows only the latest update to the density, 
where the initial state is displayed in the background. In 
each panel, the density at the final integration step is 
plotted in black. Pronounced matter-wave interference 
is observed at short times (and continues to grow in at 
longer times), as the dynamics of the superposition state 
(jl5|) is essentially that of the two interacting pieces of a 
Schrodinger cat, with all N atoms simultaneously inter- 
fering [50]. The behavior of this superposition state is 
strikingly different from that of the previous fragmented 
state. Sensitivity to the small (tt/4)/Ni orbital phase 
offset is easily seen. 

Fig. 2] presents the orbitals <^fe(x, t) and expansion 
coefficients (t) that are appropriate for the number 
entangled state (|48|) . These orbitals and expansion co- 
efficients were obtained, from first principles, with our 
time-independent symmetry-breaking MCBSCF theory 
[Io| . As discussed previously, the determination of cor- 
rect and mutually self-consistent initial orbitals and ex- 
pansion coefficients is, by no means, an easy task. Note 
that, in distinction to the fragmented state orbitals corre- 
sponding to N/2 atoms in each well, the entangled state 
orbitals reflect the fact almost all N atoms are in both 
wells simultaneously. Consequently, the entangled state 
orbitals are much broader than the fragmented state or- 
bitals, due to their mutually repulsive interactions. This 
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FIG. 2: (Color online) Orbitals and expansion coefficients corresponding to the ballistic expansion of an initially fragmented 
ground state of the double-well condensate computed within the TDMCBSCF theory. The upper three panels display the 
dynamics of an initially stationary state of the BEC, while the lower three panels display that of an initially (n /4)/(N/2) 
orbital phase offset between the left and right wells. Time increases from bottom to top in each panel, with snapshots taken 
every 0.25 ms, so that a history of the dynamics is captured at 1 ms and 2 ms after release of the trap at t = 0. Colors alternate 
to aid in visualization with red, blue, and black representing the real, imaginary and square modulus of the fragmented state 
orbitals and expansion coefficients. The units of the orbitals and expansion coefficients are arbitrary. 



broadening manifests itself in the interference patterns of 
the number entangled state. 

3. Ballistic expansion dynamics within time- dependent 
Gross-Pitaevskii theory 

To further place the dynamics of these two con- 
densate states into perspective, the single-orbital (and 
single-configurational) time-dependent Gross-Pitaevskii 
(TDGP) equation [H|, H3| 

ih(j)(x,t) — h(x)(j)(x,t) 

+ (N — 1)[J 0* (x', i)V(x, x')0(x', i)dV]0(x, t) 

(49) 

is integrated, at both zero and 7r/4 initial phase offsets, 
to mimic the t = density distribution of the fragmented 
ground state with N = 100. Note that this model is not 
an appropriate approximation; it is impossible, within 
the TDGP formalism to accurately describe a fragmented 
state, since, by definition, a fragmented state has macro- 
scopic occupation in two or more Fock states [26|, |33j, [45J . 
Since the TDGP dynamics involves only a single con- 
figuration, at all times, it would not be surprising for 



differences, perhaps even substantial differences, to arise 
between its dynamics and the dynamics associated with 
the TDMCBSCF theory, where even an initially single- 
configurational state may, in time, evolve into a com- 
plicated superposition of many configurations. Fig. O 
presents the expansion dynamics of the TDGP approxi- 
mation to the fragmented ground state density after re- 
lease of the trapping potential. The upper three pan- 
els display the dynamics of an initially stationary state 
density, while the lower three panels display that of a 
TDGP density which has an initial phase offset of 7r/4 
between the left and right halves of the TDGP orbital 
4>. Snapshots are taken, from left to right, every 0.25 ms 
so that the dynamics is captured at 1 ms, 2 ms, and 3 
ms following the trap's release. For clarity, the histories 
ending at 2 ms and 3 ms show only a few time updates 
to the density, where, at 3 ms, the initial state is dis- 
played in the background. In each panel, the density at 
the final integration step is plotted in black. Significant 
interference occurs, with near 100% visibility, within this 
single-configurational approximation: a result that was 
demonstrated in the linear regime, over a decade ago, 
in Ref. [46|]. This severe interference is in disagreement 
with the full TDMCBSCF theory presented in Fig. [TJ 
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FIG. 3: (Color online) Ballistic expansion dynamics of a macroscopic quantum superposition state [the entangled number state 
given in Eq. (|48[) ] of the double-well condensate computed within the TDMCBSCF theory. The upper three panels display the 
dynamics of an initially stationary state, while the lower three panels display that of an initially (tt/4)/N-i orbital phase offset 
between the left and right wells, where, in this case N\ = 90. Time increases from left to right, with snapshots taken every 0.25 
ms so that a history of the dynamics is captured at 1 ms, 2 ms, and 3 ms after release of the trap at t = 0. It might appear 
that only three snapshots exist in frame one, but the color version will show that the density at t — 0.25 ms overlaps the initial 
density. Colors alternate to aid in visualization. It is strikingly apparent that interference occurs in the short time dynamics 
of the superposition state density, and that the expected sensitivity to the small (n/4:)/Ni phase offset is easily observed. 



where, as in the MIT interference experiment [43j, data 
is taken in the high density regime where neither nonlin- 
earity nor correlations can be neglected and interference 
fringes are subsequently observed with much less visibil- 
ity. Here we see that atomic correlations almost quench 
the purely mean-field effects seen in TDGP theory. It is 
our opinion that this is a previously unappreciated result. 



V. CONCLUSION 

We have presented a first principles time-dependent 
quantum many-body theory of identical bosons to de- 
scribe the many-body dynamics of a double-well BEC 
at zero temperature. Within a restricted two-state Fock 
space, our TDMCBSCF theory, which is derived from 
the TDVP, includes the full effects of the condensate's 
mean field and that of atomic correlation between atoms 
in different configurations, and, additionally, enjoys the 
complete and self-consistent coupling between the expan- 
sion coefficients of each configuration and the underlying 
mean-field orbitals. The TDMCBSCF evolution equa- 
tions form a well-defined initial value problem and are an 



approximation to the exact time-dependent many-body 
Schrodinger equation. They have been implemented in 
an efficient and general numerical algorithm. In order 
to study the role of initial atomic correlation and mean- 
field effects upon the dynamics, we explore the ballistic 
expansion of an initially fragmented condensate ground 
state and an initially macroscopic quantum superposi- 
tion state (also called a Schrodinger cat state) of the 
double- well BEC, with two different initial phase offsets 
between the left and right wells, following release of the 
external trapping potential. Brilliant matter-wave inter- 
ference is observed at all times in the expansion dynamics 
of both the fragmented ground state and the superposi- 
tion state, which is analogous to the double slit with all N 
atoms interfering simultaneously. Remarkable sensitivity 
to small phase offsets between the left and right halves of 
the macroscopic quantum superposition are manifested 
in observable shifts in their interference patterns. It is 
shown that naively approximating the dynamics of the 
fragmented state with that of the TDGP approach leads 
to drastically different results than what is obtained in 
the full many-body theory and what is observed in ex- 
periment. 
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FIG. 4: (Color online) Orbitals and expansion coefficients corresponding to the ballistic expansion of a macroscopic quantum 
superposition state [the entangled number state given in Eq. (|48|) ] of the double- well condensate computed within the TDMCB- 
SCF theory. The upper three panels display the dynamics of an initially stationary state, while the lower three panels display 
that of an initially (tv/4)/Ni orbital phase offset between the left and right wells, where, in this case Ni = 90. Time increases 
from bottom to top, with snapshots taken every 0.25 ms so that a history of the dynamics is captured at 1 ms, 2 ms, and 3 ms 
(for the expansion coefficients only) after release of the trap at t = 0. Colors alternate to aid in visualization with red, blue, 
and black representing the real, imaginary and square modulus of the fragmented state orbitals and expansion coefficients. The 
units of the orbitals and expansion coefficients are arbitrary. 



We note that, of course, decoherence will need to be 
taken into account [53j |. especially for high- lying macro- 
scopic quantum superposition states [54| . In addition, if 
observations are being made then quantum back-action 
should also be included (5f| However, getting the 

many-body quantum physics correct in the absence of 
these effects is a prerequisite to the proper treatment of 
decoherence and observation. The former is the contri- 
bution of this paper. 

While in the final preparation of this paper, an article 
of similar intent has appeared on the arXiv; see Ref. [57| . 
However, the short length of Ref. [13] has prevented the 
presentation of any real details and, consequently, it is 
not possible to make a detailed comparison to the present 
work. Furthermore, we note that the numerical examples 



presented here are quite different from the fragmentation 
example in Ref. [57] . 
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